% numjac.m

% take numerical partial derivatives of the function, y=f(x,...),
% where y is a n by 1 (column) vector and x is a m x 1 vector. 

function D=numjac(fcnname,x,delta,z1,z2,z3,z4)

m=length(x);

ybar=feval(fcnname,x,z1,z2,z3,z4);

D=[];

% create the discrete differences, propoportional
% to the size of x. More specifically, create
% a matrix of f(i,j), with this being the derivative
% of the ith function with respect to the jth variable

for j=1:m;
   xnew=x;
   step=x(j)*delta;
   xnew(j)=xnew(j)+step;
   ynew=feval(fcnname,xnew,z1,z2,z3,z4);
   dy=(ynew-ybar)/step;
   D=[D dy]; 
end
